Obliczenie tempa reakcji termojądrowych wymaga zrozumienia trzech składowych rachunku:

1. Statystyczne uśrednianie w gazie
2. Tunelowanie przez barierę kulombowską
3. Obliczenie czynników w przekroju czynnym, wynikających z konkretnego oddziaływania.


Poniżej naszkicuję wyprowadzenie czynnika, znanego jako astrofizyczny S-factor, dla pierwszej, kluczowej reakcji cyklu pp (proton+proton):




Powyższa reakcja zachodzi poprzez oddzialywanie słabe, o czym świadczy obecność neutrina elektronowego w kanale wyjściowym. Ogólną teorią takich oddziaływań jest model elektroslaby, będący częścią Modelu Standardowego. Na nasze potrzeby zderzenia przy badzo niskich energiach w zupełności wystarcza jednak  ,,stary’’ model Fermiego z oddziaływaniem czteropunktowym, bez bozonów pośredniczących . W takim modelu, reakcja pp jest bliźniaczo podobna do reakcji rozpadu neutronu, która to reakcja jest jedną z ,,najprostszych’’ reakcji słabych - prostszy jest tylko rozpad mionu, który w odróżnieniu od neutronu (złożonego z kwarków) jest cząstką ,,punktową’’ (elementarną). 

Odpowiedni diagram Feynmanna jest więc ten sam, o ile potraktujemy reakcję pp jak ,,rozpad diprotonu’’. Różne będą tylko masy i ich różnice (a więc ,,rozmiar’’ przestrzeni fazowej), natomiast element macierzowy pozostaje ten sam. 



VS

Ponieważ czas rozpadu neutronu jest ,,znany’’, można go użyć jako stałej proporcojonalności, zamiast stałej Fermiego - siły sprzężenia oddziaływań słabych. ( NA podstawie THE RATE OF THE PROTON-PROTON REACTION AND SOME RELATED REACTIONS John N. Bahcall and Robert M. May The Astrophysical Journal, Vol. 155, February 1969). Zaczniemy od wyznaczenia tzw. współczynnika przestrzeni fazowej f. 

Rachunek pokazuję w całości, gdyż równocześnie jest on wyprowadzeniem przybliżonych wzorów na teoretyczny kształt widma neutrn słonecznych, a przynajmniej jego ciągłej części wynikającej z rozpadów beta. 


W wyrażeniach na tempo rozpadu pojawiają się całki po pędach wyjściowych cząstek, tj. elektronów  i neutrin :



 Przechodzimy do układu sferycznego i całkując po wszystkich kierunkach elektronów i neutrin mamy:
 
 
 
 Dla elektronów , dla neutrin  (c=1). Całkowanie delty Diraca
sprowadza się do wyrugowania energii elektronu, bądź neutrina. W teorii rozpadu $$
na ogół ruguje się , my robimy odwrotnie:




czyli całkowite tempo rozpadu jest proporcjonalne do:

Po zamianie zmiennych na bezwymiarowe     mamy:

Integrate[(q-x)x^2Sqrt[(q-x)^2-1],{x,0,q-1},Assumptions->q>1]

Plot[(1/60) ((Sqrt[-1+(q^2)]) (-8-9 (q^2)+2 (q^4))+15 q Log[q+(Sqrt[-1+(q^2)])]),{q,1,2}]

me=Quantity[1,"ElectronMass"];

Mp=Quantity[1,"ProtonMass"];

Md=Quantity[1,"DeuteronMass"];

Mn=Quantity[1,"NeutronMass"]

Obliczamy współczynniki dla neutronu i ,,diprotonu’’:

τ=Quantity[613.9,"Seconds"]

fn=(1/60) ((Sqrt[-1+(q^2)]) (-8-9 (q^2)+2 (q^4))+15 q Log[q+(Sqrt[-1+(q^2)])])/.q->(Mn-Mp)/me

fpp=(1/60) ((Sqrt[-1+(q^2)]) (-8-9 (q^2)+2 (q^4))+15 q Log[q+(Sqrt[-1+(q^2)])])/.q->(2Mp-Md)/me

c=Quantity[1,"SpeedOfLight"];

α=Quantity[1,"FineStructureConstant"]

Poniżej zapisano podręcznikowe wyniki dla rozpadu neutronu, ignorujące detale struktury nukleonów.

rate1=16 Pi^2 2GF0^2*(me c^2)^5fn/(2Pi)^6*2Pi/hbar*2

rate2=2 fn/Pi^3/hbar^7*(GF0*hbar^3 c^3)^2 me^5 c^4 (*Neutron Beta-Decay Christopher B.Hayes December 6,2012 *)

rate1/rate2

Dosyć powszechne (cyt. Bahcalla: “It is both conventional and convenient”, The Astrophysical Journal, Vol. 155, February 1969) jest zastąpienie przekroju czynnego astrofizycznym S-factorem:

 


Wynik referencyjny, do którego będziemy ,,dążyć’’ to

S0=Quantity[4.07*10^-22 ,"Kiloelectronvolts"*"Barns"] (* Zmierzony S-factor *)

Λ=Sqrt[7]

(Entity["deuteron", Entity["Particle","Deuteron"], "Entity[\"Particle\", \"Deuteron\"]", "particle"])[(EntityProperty["classical radius", EntityProperty["Particle","ClassicalRadius"], "EntityProperty[\"Particle\", \"ClassicalRadius\"]"])]

r0=Quantity[2,"Femtometers"]

hbar=Quantity[1,"ReducedPlanckConstant"]

Sqrt[Mp c^2* Quantity[2.2,"MeV"]]/hbar/c//UnitConvert

1/%

4 Pi^2 c α Mp Log[2]*fpp/(fn*τ) * Λ^2*(2r0)^3

UnitConvert[%]

SBahcall=UnitConvert[%,"Barns*MeV"]

S0/SBahcall

fn τ

fpp